#include "f2c.h"
doublereal duni_0_(int n__, integer *iseed);
doublereal duni_(void);
doublereal dustar_(integer *iseed);
doublereal formf_(doublereal *p);
/* Subroutine */ int matin1_(doublereal *a, integer *dim1, integer *n1, 
	integer *n2, integer *indexx, integer *nerror, doublereal *determ)
{
    /* System generated locals */
    integer i__1, i__2, i__3;
    doublereal d__1;

    /* Local variables */
    static integer i__, n, i1, i2, i3, dim, main, emat, icol, jcol, lpiv;
    static doublereal swap;
    static integer pcol1, pcol2, nmin1;
    static doublereal deter, pivot;
    static integer pivcol;

/* ---- LINEAR SYSTEM SOLVER FROM CERN LIBRARY. */
/* ---- CALLED BY EMEAS, NLS */
    /* Parameter adjustments */
    --indexx;
    --a;

    /* Function Body */
    deter = 1.f;
    n = *n1;
    emat = n + *n2;
    dim = *dim1;
    nmin1 = n - 1;
    pivcol = 1 - dim;
    i__1 = n;
    for (main = 1; main <= i__1; ++main) {
	pivot = 0.f;
	pivcol += dim;
	pcol1 = pivcol + main - 1;
	pcol2 = pivcol + nmin1;
	i__2 = pcol2;
	for (i__ = pcol1; i__ <= i__2; ++i__) {
/*     第一次循环1~N-1循环，找到A对角线绝对值最大的元素，赋值给pivot，对应的行标LPIV */
	    if ((d__1 = a[i__], abs(d__1)) - abs(pivot) <= 0.) {
		goto L20;
	    } else {
		goto L10;
	    }
L10:
	    pivot = a[i__];
	    lpiv = i__;
L20:
	    ;
	}
	if (pivot != 0.) {
	    goto L30;
	} else {
	    goto L150;
	}
/*     如果对角线上元素都是0，nerror=-1,说明det|A|太小了，stop */
L30:
	icol = lpiv - pivcol + 1;
	indexx[main] = icol;
	if (icol - main <= 0) {
	    goto L60;
	} else {
	    goto L40;
	}
L40:
	deter = -deter;
	icol -= dim;
	i3 = main - dim;
	i__2 = emat;
	for (i__ = 1; i__ <= i__2; ++i__) {
	    icol += dim;
	    i3 += dim;
	    swap = a[i3];
	    a[i3] = a[icol];
/* L50: */
	    a[icol] = swap;
	}
L60:
	deter *= pivot;
	pivot = 1.f / pivot;
	i3 = pivcol + nmin1;
	i__2 = i3;
	for (i__ = pivcol; i__ <= i__2; ++i__) {
/* L70: */
	    a[i__] = -a[i__] * pivot;
	}
	a[pcol1] = pivot;
	i1 = main - dim;
	icol = 1 - dim;
	i__2 = emat;
	for (i__ = 1; i__ <= i__2; ++i__) {
	    icol += dim;
	    i1 += dim;
	    if (i__ - main != 0) {
		goto L80;
	    } else {
		goto L100;
	    }
L80:
	    jcol = icol + nmin1;
	    swap = a[i1];
	    i3 = pivcol - 1;
	    i__3 = jcol;
	    for (i2 = icol; i2 <= i__3; ++i2) {
		++i3;
/* L90: */
		a[i2] += swap * a[i3];
	    }
	    a[i1] = swap * pivot;
L100:
	    ;
	}
/* L110: */
    }
    i__1 = n;
    for (i1 = 1; i1 <= i__1; ++i1) {
	main = n + 1 - i1;
	lpiv = indexx[main];
	if (lpiv - main != 0) {
	    goto L120;
	} else {
	    goto L140;
	}
L120:
	icol = (lpiv - 1) * dim + 1;
	jcol = icol + nmin1;
	pivcol = (main - 1) * dim + 1 - icol;
	i__2 = jcol;
	for (i2 = icol; i2 <= i__2; ++i2) {
	    i3 = i2 + pivcol;
	    swap = a[i2];
	    a[i2] = a[i3];
/* L130: */
	    a[i3] = swap;
	}
L140:
	;
    }
    *determ = deter;
    *nerror = 0;
    return 0;
L150:
    *nerror = -1;
    *determ = deter;
    return 0;
} /* matin1_ */

doublereal duni_0_(int n__, integer *iseed)
{
    /* Initialized data */

    static doublereal u[17] = { .471960981577884755837789724978,
	    .930323453205669578433639632431,.110161790933730836587127944899,
	    .057150199627313951836263875701,.402467554779738266237538503137,
	    .451181953427459489458279456915,.296076152342721102174129954053,
	    .0128202189325888116466879622359,.314274693850973603980853259266,
	    .00335521366752294932468163594171,.48868504520043937160785036784,
	    .195470426865656758693860613516,.864162706791773556901599326053,
	    .335505955815259203596381170316,.377190200199058085469526470541,
	    .400780392114818314671676525916,.374224214182207466262750307281 };
    static integer i__ = 17;
    static integer j = 5;
    static integer k = 47;
    static doublereal c__ = .065104166666593244;

    /* System generated locals */
    integer i__1;
    doublereal ret_val;

    /* Local variables */
    static doublereal s, t;
    static integer i1, j1, k1, l1, m1, ii, jj;

/* ***BEGIN PROLOGUE  DUNI */
/* ***DATE WRITTEN   880714 (YYMMDD) */
/* ***REVISION DATE  880714 (YYMMDD) */
/* ***KEYWORDS  RANDOM NUMBERS, UNIFORM RANDOM NUMBERS */
/* ***AUTHOR    KAHANER, DAVID, SCIENTIFIC COMPUTING DIVISION, NBS */
/*             MARSAGLIA, GEORGE, SUPERCOMPUTER RES. INST., FLORIDA ST. U. */

/* ***PURPOSE  THIS ROUTINE GENERATES DOUBLE PRECISION UNIFORM */
/*             RANDOM NUMBERS ON [0,1) */
/* ***DESCRIPTION */
/*        COMPUTES DOUBLE PRECISION UNIFORM NUMBERS ON [0,1). */
/*           FROM THE BOOK, "NUMERICAL METHODS AND SOFTWARE" BY */
/*                D. KAHANER, C. MOLER, S. NASH */
/*                PRENTICE HALL, 1988 */

/*       USAGE: */
/*              TO INITIALIZE THE GENERATOR */
/*                   USEED = DUSTAR(ISEED) */
/*               WHERE: ISEED IS ANY NONZERO INTEGER */
/*                  WILL RETURN FLOATING POINT VALUE OF ISEED. */

/*               SUBSEQUENTLY */
/*                       U = DUNI() */
/*                  WILL RETURN A REAL UNIFORM ON [0,1) */

/*                ONE INITIALIZATION IS NECESSARY, BUT ANY NUMBER OF EVALUATIONS */
/*                  OF DUNI IN ANY ORDER, ARE ALLOWED. */

/*           NOTE: DEPENDING UPON THE VALUE OF K (SEE BELOW), THE OUTPUT */
/*                       OF DUNI MAY DIFFER FROM ONE MACHINE TO ANOTHER. */

/*           TYPICAL USAGE: */

/*               DOUBLE PRECISION U,DUNI,DUSTAR,USEED */
/*               INTEGER ISEED */
/* C                 SET SEED */
/*               ISEED = 305 */
/*               USEED = DUSTAR(ISEED) */
/*               DO 1 I = 1,1000 */
/*                   U = DUNI() */
/*             1 CONTINUE */
/* C                 NOTE: IF K=47 (THE DEFAULT, SEE BELOW) THE OUTPUT VALUE OF */
/* C                           U WILL BE 0.812053811384E-01... */
/*               WRITE(*,*) U */
/*               END */


/* ***REFERENCES  MARSAGLIA G., "COMMENTS ON THE PERFECT UNIFORM RANDOM */
/*                 NUMBER GENERATOR", UNPUBLISHED NOTES, WASH S. U. */
/* ***ROUTINES CALLED  (NONE) */
/* ***END PROLOGUE DUNI */
/*                            2**47=0.140737488355328D+15 */

/*      LOAD DATA ARRAY IN CASE USER FORGETS TO INITIALIZE. */
/*      THIS ARRAY IS THE RESULT OF CALLING DUNI 100000 TIMES */
/*         WITH ISEED=305 AND K=96. */
    switch(n__) {
	case 1: goto L_dustar;
	}


/*   BASIC GENERATOR IS FIBONACCI */

    ret_val = u[(0 + (0 + (i__ - 1 << 3))) / 8] - u[(0 + (0 + (j - 1 << 3))) /
	     8];
    if (ret_val < 0.) {
	ret_val += 1.;
    }
    u[i__ - 1] = ret_val;
    --i__;
    if (i__ == 0) {
	i__ = 17;
    }
    --j;
    if (j == 0) {
	j = 17;
    }

/*   SECOND GENERATOR IS CONGRUENTIAL */

    c__ += -.54387223503964321;
    if (c__ < 0.) {
	c__ += .99999999999918288;
    }

/*   COMBINATION GENERATOR */

    ret_val -= c__;
    if (ret_val < 0.) {
	ret_val += 1.;
    }
    return ret_val;


L_dustar:

/*          SET UP ... */
/*          CONVERT ISEED TO FOUR SMALLISH POSITIVE INTEGERS. */

    i1 = abs(*iseed) % 177 + 1;
    j1 = abs(*iseed) % 167 + 1;
    k1 = abs(*iseed) % 157 + 1;
    l1 = abs(*iseed) % 147 + 1;

/*              GENERATE RANDOM BIT PATTERN IN ARRAY BASED ON GIVEN SEED. */

    for (ii = 1; ii <= 17; ++ii) {
	s = 0.;
	t = .5;
/*             DO FOR EACH OF THE BITS OF MANTISSA OF WORD */
/*             LOOP  OVER K BITS, WHERE K IS DEFAULTED TO 47 */
	i__1 = k;
	for (jj = 1; jj <= i__1; ++jj) {
	    m1 = i1 * j1 % 179 * k1 % 179;
	    i1 = j1;
	    j1 = k1;
	    k1 = m1;
	    l1 = (l1 * 53 + 1) % 169;
	    if (l1 * m1 % 64 >= 32) {
		s += t;
	    }
/* L3: */
	    t *= .5;
	}
/* L2: */
	u[ii - 1] = s;
    }
    ret_val = (real) (*iseed);
    return ret_val;
} /* duni_ */

doublereal duni_(void)
{
    return duni_0_(0, (integer *)0);
    }

doublereal dustar_(integer *iseed)
{
    return duni_0_(1, iseed);
    }
doublereal formf_(doublereal *p)
{
    /* Initialized data */

    static doublereal ftbl[42]	/* was [21][2] */ = { 1.,.926,.861,.803,.75,
	    .704,.661,.623,.588,.556,.527,.5,.476,.453,.432,.413,.394,.378,
	    .362,.347,.333,0.,.0067,.0203,.0372,.0558,.0754,.0954,.115,.135,
	    .155,.174,.192,.21,.227,.244,.26,.276,.291,.306,.32,.333 };

    /* System generated locals */
    doublereal ret_val, d__1;

    /* Builtin functions */
    double sqrt(doublereal), log(doublereal);

    /* Local variables */
    static doublereal a, d__;
    static integer i__, j;
    static doublereal q, dsq;

/* ---- CALCULATE THE SPACE-CHARGE FORM FACTOR FOR ELLIPSOIDAL BEAMS */
/* ---- CALLED BY SPCHG */
    q = abs(*p);
    j = 1;
    if (q <= 1.f) {
	goto L5;
    }
    if (q > 20.f) {
	goto L10;
    }
    j = 2;
    q = 1.f / q;
L5:
    i__ = (integer) (q * 20.f + 1.f);
    a = q * 20.f - (real) (i__ - 1);
    ret_val = (1.f - a) * ftbl[i__ + j * 21 - 22] + a * ftbl[i__ + 1 + j * 21 
	    - 22];
    return ret_val;
L10:
/* Computing 2nd power */
    d__1 = q;
    dsq = d__1 * d__1 - 1.f;
    d__ = sqrt(dsq);
    ret_val = (q * log((real) (q + d__)) / d__ - 1.f) / dsq;
    return ret_val;
} /* formf_ */

